clear
clc
close all

load RTDB_NAIRU

time=datevec(dates);
start=find(time(:,1)==1995 & time(:,2)==1);
eend=find(time(:,1)==2009 & time(:,2)==4);

% transform data
HICP       =  HICP(start:eend,:);%250
UR         =  UR(start:eend,:);%313
OIL        =  OIL(start:eend,:);%193
dates      =  dates(start:eend,:);

inf    =  (log(HICP(13:end,:))-log(HICP(1:end-12,:)))*100;
oil    =  (log(OIL(13:end,:))-log(OIL(1:end-12,:)))*100;
UR     = UR(13:end,:);
dates  =  dates(13:end,:);

dinf   =  inf(2:end,:)-inf(1:end-1,:);
doil   =  oil(2:end,:)-oil(1:end-1,:);
UR     =  UR(2:end,:);
dates  =  dates(2:end,:);

UR_2   = UR(1:end-2,:);
UR_1   = UR(2:end-1,:);
UR   = UR(3:end,:);

dinf_1 =  dinf(2:end-1,:);
dinf   =  dinf(3:end,:);

doil_1 =  doil(2:end-1,:);
doil   =  doil(3:end,:);

dates  =  dates(3:end,:);
time=datevec(dates);

[T,nv]=size(doil);

NAIRU=[];Band=[];NAIRU_pseudo=[];Test=[];U=[];TR=[];

% compute NAIRU
for i=1:nv    
    reg      =  [ones(T,1) UR_1(:,i) UR_2(:,i)...
            dinf_1(:,i)  doil_1(:,i) ];
    mm=min(sum(isfinite([reg dinf(:,i)])));
    reg=reg(1:mm,:);
    
    reg_pseudo=[ones(mm,1) UR_1(1:mm,end) UR_2(1:mm,end)...
            dinf_1(1:mm,end)  doil_1(1:mm,end) ];
    
    coef     =  inv(reg' *reg)*reg'*dinf(1:mm,i);
    
    coef_pseudo=  inv(reg_pseudo' *reg_pseudo)*reg_pseudo'*dinf(1:mm,end);
    
    NAIRU    =  [NAIRU; -coef(1)/(coef(2)+coef(3))];
    
    NAIRU_pseudo    =  [NAIRU_pseudo; -coef_pseudo(1)/(coef_pseudo(2)+coef_pseudo(3))];
    
    for u = 5:.1:15
        reg= [ones(mm,1) UR_1(1:mm,i)-u UR_2(1:mm,i)-u dinf_1(1:mm,i) doil(1:mm,i) ];
        coef     =  inv(reg' *reg)*reg'*dinf(1:mm,i);
        e        =  dinf(1:mm,i)-reg*coef;
        v        =  var(e);
        S        =  inv(reg' *reg);
        tr       =  coef(1)/sqrt(v*S(1,1));
        if abs(tr)<=1.645
            Test=[Test 1];
        else
            Test=[Test 0];
        end
        U=[U u];
        TR=[TR tr];
    end
    Band=[Band;min(U(Test==1)) max(U(Test==1))];
    Test=[];U=[];TR=[];
end

% plot NAIRU in real time and pseudo real time
figure(1)
plot(datenum(VV),Band(:,1),'-.r',datenum(VV),NAIRU,'*-b',datenum(VV),NAIRU_pseudo,'d-g',datenum(VV),ones(nv,1)*NAIRU(end),'--m',datenum(VV),Band(:,2),'-.r','LineWidth',2),...
    datetick('x',12),legend('90% c.i.','RT NAIRU','Pseudo NAIRU','Last NAIRU'),axis([datenum(VV(1,:)) datenum(VV(end,:))  4 12]),xlabel('vintages')


%compute unemployment gap
T=length(UR);
UR_gap=UR-ones(T,1)*NAIRU';
last_obs=nan(T,1);

lo=find(isfinite(UR_gap(:,1)));

for i=lo(end):T
    sl=find(isfinite(UR_gap(i,:)));
    last_obs(i,:)=UR_gap(i,sl(1));
end

% plot unemployment gap

figure(2)
plot(dates(1:157,:),UR_gap(1:157,:),dates(1:157,:),last_obs,'x-r','LineWidth',2,'MarkerSize',15),...
    datetick('x',12),xlabel('time')

